clear
rng default
n_types=4;
all_years=1940:1:1967;
n_years=size(all_years,2);

%% Construct pictures reporting the bidimensional projections of the identified sets computed in the cluster. 
%When the hitting the grid bounds report -Inf, +Inf in the table below. 
%3 indicates -Inf, +Inf
%2 indicates lower bound -Inf
%1 indicates upper bound +Inf
%NaN indicates finite lower and upper bound

%U	2	3	4	2	3	4	2	3	4	2	3	4
U=[ 1	3	3	1	3	3	1	3	3	1	1	1; ...
 	1	3	3	1	3	3	1	3	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	1];

%V	2	3	4	2	3	4	2	3	4	2	3	4
V=[ NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1];


U_reshape=cell(n_years,1);
for h=1:n_years
    U_reshape{h}=NaN(n_types, n_types);
    for i=1:n_types
        U_reshape{h}(i,2:n_types)=U(h,(i-1)*(n_types-1)+1:i*(n_types-1));
    end
end

V_reshape=cell(n_years,1);
for h=1:n_years
    V_reshape{h}=NaN(n_types, n_types);
    for i=1:n_types
        V_reshape{h}(2:n_types,i)=(V(h,(i-1)*(n_types-1)+1:i*(n_types-1))).';
    end
end

%% Logit estimates

Phi_Gumbel_weighted=cell(n_years,1);
weight=cell(n_years,1);
PYX_cond=cell(n_years,1);
PXY_cond=cell(n_years,1);

for h=1:n_years
    
year=all_years(h);

% Weights
weight_temp=load(['weight' num2str(h) '_4g']);
weight{h}=weight_temp.x;

% Probability distribution of Ychosen|X
condprobsmen_temp=load('probs_allyears_data.mat',['condprobsmen' num2str(year) '_140']);
condprobsmen=condprobsmen_temp.(['condprobsmen' num2str(year) '_140']);
PYX_cond{h}=[condprobsmen(:,2:end) condprobsmen(:,1)]; %4x5
%Note: PYX_cond(x,y)=Probability to marry woman y born in any year conditional on being man x born in 1940

% Probability distribution of Xchosen|Y
condprobswomen_temp=load('probs_allyears_data.mat',['condprobswomen' num2str(year+1) '_140']);
condprobswomen=condprobswomen_temp.(['condprobswomen' num2str(year+1) '_140']);
PXY_cond{h}=[condprobswomen(:,2:end) condprobswomen(:,1)];  %4x5
%Note: PXY_cond(y,x)=Probability to marry man x born in any year conditional on being woman y born in 1941

Phi_Gumbel_temp=param_Gumbel_closedform_nosingles(PYX_cond{h}, PXY_cond{h}, n_types, n_types);

Phi_Gumbel_weighted{h}=Phi_Gumbel_temp.*weight{h};

end

Phi_Gumbel_early=(Phi_Gumbel_weighted{1}+Phi_Gumbel_weighted{2}+Phi_Gumbel_weighted{3})./(weight{1}+weight{2}+weight{3});

Phi_Gumbel_late=(Phi_Gumbel_weighted{26}+Phi_Gumbel_weighted{27}+Phi_Gumbel_weighted{28})./(weight{26}+weight{27}+weight{28});

Delta_Gumbel_Phi=Phi_Gumbel_late-Phi_Gumbel_early;

%% Our estimates

fp='.../SpecB/output'; %specify path
Phi_weighted=cell(n_years,1); 


for h=1:n_years
    Phi_weighted{h}=cell(n_types, n_types);
    
    IdSetM=cell(n_types,1);
    IdSetW=cell(n_types,1);
    
    for j=1:n_types
        
       filepath = fullfile( fp, ['IdSetM' num2str(all_years(h)) '_final_' num2str(j) '.mat'] );
       if exist( filepath, 'file' )
          A = getfield( load( filepath ), 'IdSetM_final' ); 
          IdSetM{j}=A;
       end 
       
       filepath = fullfile( fp, ['IdSetW' num2str(all_years(h)) '_final_' num2str(j) '.mat'] );
       if exist( filepath, 'file' )
          A = getfield( load( filepath ), 'IdSetW_final' ); 
          IdSetW{j}=A;
       end  
    end
    
    
for i=2:n_types+1
    
for j=2:n_types+1

    M=unique(IdSetM{i-1}(:,j), 'rows', 'stable');  
    W=unique(IdSetW{j-1}(:,i), 'rows', 'stable');  
    
    sm=size(M,1);
    temp=cell(sm,1);
    temp_V=cell(sm,1);
    for k=1:sm
        temp{k}=M(k)+W(:); 
    end
    Phi_temp=vertcat(temp{:});
    
    Phi_int=[min(Phi_temp); max(Phi_temp)];

    if  U_reshape{h}(i-1,j-1)==1 || U_reshape{h}(i-1,j-1)==3 || V_reshape{h}(i-1,j-1)==1 || V_reshape{h}(i-1,j-1)==3 
        Phi_int(2)=Inf;
    end
    if  U_reshape{h}(i-1,j-1)==2 || U_reshape{h}(i-1,j-1)==3 || V_reshape{h}(i-1,j-1)==2 || V_reshape{h}(i-1,j-1)==3 
        Phi_int(1)=-Inf;
    end
    
    Phi_int(1)=(Phi_int(1)+log(PYX_cond{h}(i-1,end))+log(PXY_cond{h}(j-1,end)))*weight{h}(i-1,j-1);
    Phi_int(2)=(Phi_int(2)+log(PYX_cond{h}(i-1,end))+log(PXY_cond{h}(j-1,end)))*weight{h}(i-1,j-1);
    
    Phi_weighted{h}{i-1,j-1}=[Phi_int(1) Phi_int(2)];
end
end
end

Phi_early_R2=cell(n_types, n_types);
Phi_late_R2=cell(n_types, n_types);
Delta_Phi_R2=cell(n_types, n_types);
for i=1:n_types
    for j=1:n_types
    Phi_early_R2{i,j}=(Phi_weighted{1}{i,j}+Phi_weighted{2}{i,j}+Phi_weighted{3}{i,j})./(weight{1}(i,j)+weight{2}(i,j)+weight{3}(i,j));

    Phi_late_R2{i,j}=(Phi_weighted{26}{i,j}+Phi_weighted{27}{i,j}+Phi_weighted{28}{i,j})./(weight{26}(i,j)+weight{27}(i,j)+weight{28}(i,j));

    Delta_Phi_R2{i,j}(1)=Phi_late_R2{i,j}(1)-Phi_early_R2{i,j}(2);
    Delta_Phi_R2{i,j}(2)=Phi_late_R2{i,j}(2)-Phi_early_R2{i,j}(1);
    end
end

save('Table1_R2', 'Phi_Gumbel_early', 'Phi_Gumbel_late', 'Delta_Gumbel_Phi', 'Phi_early_R2', 'Phi_late_R2', 'Delta_Phi_R2')
